function [cur,wrapping_size] = curvelet_limitedAngle_outerFading(imSize,radial,angular, lim_angl)
%custom curvelet-like implementation adapted to a limited angle problem.
%Fading happens autside the given angle.
%
%author: Sebastian Schmelcher; version: alpha

identifier=['limitedAngle_outerFading_' num2str(lim_angl)];
map=get_indexStruct(identifier,imSize);

num_radial=length(map);

n=imSize(1);
m=imSize(2);

origin=get_origin([n,m]);




if(radial==0)
    cur=[0.928241517645832,1,0.928241517645832;1,1,1;0.928241517645832,1,0.928241517645832];
    point=[origin(1)-1 origin(2)-1; origin(1)+1 origin(2)+1 ];
    cur={cur,point};
    wrapping_size=[3, 3];
    return;
end

if(radial>num_radial)
    return;
end

if(angular>map(radial))
    return;
end



curvelet_width=2*lim_angl/(map(radial)-2/3);
lim_angl=lim_angl+1/3*curvelet_width;
curvelet_dir=-lim_angl/2+(angular-.5)*curvelet_width;
if(angular>map(radial)/2)
    curvelet_dir=curvelet_dir+180-lim_angl;
end

curvelet_dir=angle(exp(1i*pi*curvelet_dir/180));
curvelet_width=pi*curvelet_width/180;
curvelet_supp_width=2/3*curvelet_width; %half of support width
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% length and width of mother curvelet %%%
l=ceil(5*2^(radial)/3)-floor(2^(radial+1)/3*cos(curvelet_supp_width));
b=2*ceil(abs(5/3*2^(radial)*sin(curvelet_supp_width)));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% compute rectangle containg complete curvelet %%%   
r1=floor(2^(radial+1)/3);%*cos(curvelet_supp_width));
r2=ceil(5*2^(radial)/3);
if(abs(curvelet_dir) < 3/4*pi&& abs(curvelet_dir) >pi/4)
    %North or South
   
    x=r2*[cos(curvelet_dir+curvelet_supp_width),cos(curvelet_dir-curvelet_supp_width)];        
    x=[x,r1*[cos(curvelet_dir+curvelet_supp_width),cos(curvelet_dir-curvelet_supp_width)]];      
    y=[r1*[sin(curvelet_dir+curvelet_supp_width),sin(curvelet_dir-curvelet_supp_width)],r2*[sin(curvelet_dir+curvelet_supp_width),sin(curvelet_dir-curvelet_supp_width)]];
    if( abs(curvelet_dir-pi/2)<curvelet_supp_width/2 || abs(curvelet_dir-pi/2)>pi-curvelet_supp_width/2)
        %curvelet intersects vertical axis
        if(curvelet_dir<0)
            y(end+1)=-r2;
        else
            y(end+1)=r2; 
        end
    end
else
    %East or West  

    x=r2*[cos(curvelet_dir+curvelet_supp_width),cos(curvelet_dir-curvelet_supp_width)];        
    x=[x,r1*[cos(curvelet_dir+curvelet_supp_width),cos(curvelet_dir-curvelet_supp_width)]];    
    y=[r1*[sin(curvelet_dir+curvelet_supp_width),sin(curvelet_dir-curvelet_supp_width)],r2*[sin(curvelet_dir+curvelet_supp_width),sin(curvelet_dir-curvelet_supp_width)]];
    if( abs(curvelet_dir)<curvelet_supp_width/2 || abs(curvelet_dir)>pi-curvelet_supp_width/2)
        %curvelet intersects horizontal axis
        if( abs(curvelet_dir)<pi/2)
            x(end+1)=r2;
        else
            x(end+1)=-r2;
        end
    end
end
point_ul=[max(1,origin(1)+floor(min(y))),max(1,origin(2)+floor(min(x)))];
point_lr=[min(n,origin(1)+ceil(max(y))),min(m,origin(2)+ceil(max(x)))];
point=[point_ul;point_lr];



%%%%%%%%%%% original %%%%%%%%%%%
% [X,Y]=meshgrid(1-origin(1):n-origin(1),1-origin(2):m-origin(2));
% [theta,rad]=cart2pol(X,Y);
% rad=(2^(-radial)).*rad;
% theta=angle(exp(1i*(theta-curvelet_dir)))/curvelet_width*2;
% if(radial==num_radial)
%     cur1=lastRadialWindow(rad).*angularWindow(theta); %no upper limit
% else
%     cur1=radialWindow(rad).*angularWindow(theta);
% end

%%%%%%%%% faster %%%%%%%%%%%%
[X,Y]=meshgrid(point(3)-origin(1):point(4)-origin(1),point(1)-origin(2):point(2)-origin(2));
[theta,rad]=cart2pol(X,Y);
rad=(2^(-radial)).*rad;
theta=angle(exp(1i*(theta-curvelet_dir)))/curvelet_width*2;
if(radial==num_radial)
    cur2=lastRadialWindow(rad).*angularWindow(theta); %no upper limit
else
    cur2=radialWindow(rad).*angularWindow(theta);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% compute size of wrapping rectangle %%%
a=ceil(abs(b*sin(curvelet_dir)));
%curvelet_dir
if(abs(curvelet_dir) < 3/4*pi && abs(curvelet_dir) >pi/4)
    %North or South
    wrapping_size=[ceil(abs(l*sin(curvelet_dir)))+a,ceil(abs(b/sin(curvelet_dir)))];
else
    %East or West
    wrapping_size=[ceil(abs(b/cos(curvelet_dir))),ceil(abs(l*cos(curvelet_dir)))+a];
end
%anz_coef=prod(wrapping_size)
%cur(point(1):point(2),point(3):point(4))=cur2;
cur={cur2,point};



% -------------------------------------------------------------------------
function x = lastRadialWindow(r)
%no upper limit
rSize=size(r);

x= (r>=(5/6));
x=x+cos((5*ones(rSize)-6*r)*(pi/2)).*(r>=(2/3)) .* (r<(5/6));







% -------------------------------------------------------------------------

function x = angularWindow(t)
t=abs(t);
tSize=size(t);


x=(t<=2/3);
x=x+((t<=4/3).*(t>2/3) .* cos(  (3/2*t-ones(tSize))*(pi/2) ));

% x=(t<=(1/3));
% x=x+((t<=(2/3)) .* (t>1/3)).*cos(  (3*t-ones(tSize))*(pi/2) );



% -------------------------------------------------------------------------
function x = radialWindow(r)
rSize=size(r);

x=( (r>=(5/6)) .* (r<=(4/3)) );
x=x+cos((5*ones(rSize)-6*r)*(pi/2)).*(r>=(2/3)) .* (r<(5/6));
x=x+cos((3*r-4*ones(rSize))*(pi/2)).*(r>(4/3)) .* (r<=(5/3));

